Reconstruction of photon statistics using low performance photon counters 
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The output of a photodetector consists of a current pulse whose charge has the statistical distri- 
bution of the actual photon numbers convolved with a Bernoulli distribution. Photodetectors are 
characterized by a nonunit quantum efficiency, i.e. not all the photons lead to a charge, and by a 
finite resolution, i.e. a different number of detected photons leads to a discriminable values of the 
charge only up to a maximum value. We present a detailed comparison, based on Monte Carlo sim- 
ulated experiments and real data, among the performances of detectors with different upper limits of 
counting capability. In our scheme the inversion of Bernoulli convolution is performed by maximum- 
likelihood methods assisted by measurements taken at different quantum efficiencies. We show that 
detectors that are only able to discriminate between zero, one and more than one detected photons 
are generally enough to provide a reliable reconstruction of the photon statistics for single-peaked 
distributions, while detectors with higher resolution limits do not lead to further improvements. In 
addition, we demonstrate that, for semiclassical states, even on/off detectors are enough to provide a 
good reconstruction. Finally, we show that a reliable reconstruction of multi-peaked distributions re- 
quires either higher quantum efficiency or better capability in discriminating high number of detected 
photons. 

I. INTRODUCTION 

Reconstruction of the photon statistics, g n , of optical states provides fundamental information on the na- 
ture of any optical field and finds relevant applications in foundations of quantum mechanics, quantum state 
engineering by postselection fl], quantum information 0], and quantum metrology. Indeed, detectors with 
the capability of counting photons 0, 01 are currently under investigation. Among these, photomultiplier 
tubes (PMT's) Q] and hybrid photodetectors 0,0] are promising devices, though they have the drawback 
of a low quantum efficiency. On the other hand, solid state detectors with internal gain are still under devel- 
opment. Highly efficient thermal detectors have also been used as photon counters, though their operating 
conditions are still extreme to allow common use 0,0]. The advent of quantum tomography provided an 
alternative method to measure photon number distributions 0]. However, tomography needs the implemen- 
tation of homodyne detection, which involves challenging mode matching, especially in the case of pulsed 
optical fields. 

In principle, in a photodetector each photon ionizes a single atom, and the resulting charge is amplified 
to produce a measurable pulse. In practice, however, available photodetectors are usually characterized 
by a quantum efficiency lower than unity, which means that only a fraction of the incoming photons leads 
to an electric pulse. If the resulting current is proportional to the incoming photon flux we have a linear 
photodetector. This is, for example, the case of the high-flux photodetectors used in homodyne detection. 
On the other hand photodetectors operating at very low intensities usually resort to avalanche process or 
very high amplification in order to transform a single ionization event into a recordable pulse. Due to the 
gain instability in the process of amplification it is generally difficult to discriminate the number of detected 
photons as far as this number becomes larger. 

As a consequence the output of a photodetector consists of a current pulses whose charge statistics 
is a generalized Bernoulli convolution of the actual photon distribution with the possibility to perform a 
discrimination between the number of detected photons only up to a finite number. We define M as the 
maximum number m of detected photons that can be distinguished from m — 1, M represent the counting 
capability of the detectors: i.e. detectors that works in Geiger mode have M = 1 because form the output 
is possible to discriminate from m = 1, one detected photon, to m = 0, dark, but is not possible to 



discriminate the output for m > 2 from that for to m = 1. The inversion of such a convolution may be 
performed in several ways, though, in general it is possible only when the quantum efficiency is larger than 
77 = 0.5 iflfill . and it is inherently inefficient, as it requires a large data sample to provide a reliable result. 
On the other hand, maximum-likelihood methods assisted by measurements taken at different quantum 
efficiencies have been proved to be both effective and statistically reliable 131I1411 . In particular, it 

is possible to obtain a method to reconstruct g n without any a priori information on the state of light under 
investigation. An alternative approach is based on the so called photon chopping lll^l . i.e. on a network of 
beam splitter followed by an array of on/off detectors. Photon chopping allows to send at most one photon 
on each detector, however, with the drawback of increasing the overall complexity of the detection scheme. 

In this paper, we extend previous analysis on maximum-likelihood on/off -based schemes 13], and 
present a detailed comparison among the performances of detectors characterized by realistic quantum effi- 
ciencies and different values of M. The comparison is made using real data and by an extensive set of Monte 
Carlo simulations, performed on different states of the radiation field including quantum and semiclassical 
states associated to single-peaked as well as multi-peaked photon distributions. The reconstruction is ob- 
tained by maximum-likelihood methods assisted by measurements taken at different quantum efficiencies. 
In particular, as the statistics of the detected photons is linearly dependent on the photon distribution (see 
below), the inversion can be well approximated by using the Expectation-Maximization algorithm (EM) 
which leads to an effective and reliable iterative solution. 

Our results indicates that detectors with M = 2 discriminating between zero, one and more than one 
photon, are generally enough to provide a reliable reconstruction of the photon statistics for single-peaked 
distributions, while detectors with higher value of M do not lead to further improvements. On the other 
hand, multi-peaked distributions requires a higher quantum efficiency whereas for semiclassical states even 
on/off detectors are enough to provide a good reconstruction. 

The paper is structured as follows. In the next Section we describe the statistics of the detected photons 
and illustrate the reconstruction algorithm. In Section |ffl] we report the results of a set of Monte Carlo sim- 
ulated experiments performed on different kinds of signals, whereas in Section |W] we report experimental 
results for coherent signals. Section|V]closes the paper with some concluding remarks. 

II. STATISTICS OF DETECTED PHOTONS AND RECONSTRUCTION ALGORITHM 

Using a photodetector with quantum efficiency ?] and M = 00, that corresponding to an unlimited 
counting capability, the probability of obtaining m detected photons at the output is given by the convolution 



p v (m) = ^2 A u (m,n)g n , (1) 

n=m 

where g n = (n\g\n) is the actual photon statistics of the signal under investigation and 

A,(m,n) = ^yi-r ] r- m r ] m . (2) 

If M have a finite value only M + 1 outcomes are possible, which occur with probabilities 

q™= Pr] (m) m = 0,...,M-l (3) 

00 M-l 

< = E = 1 - E € ■ w 

m=M m=0 

Once the value of the quantum efficiency is known, Eqs. © and © provide M relations among the statistics 
of detected photons and the actual statistics of photons. At first sight the statistics of a detector with low 



counting capability (M in the range M = 1, 6) appears to provide quite a scarce piece of information 
about the state under investigation. However, if the statistics about {q™} is collected for a suitably large set 
of efficiency values, then the information is enough to reconstruct the whole photon distribution g n of the 
signal, upon a suitable truncation at N of the Hilbert space. We adopt the following strategy: by placing in 
front of the detector K filters with different transmissions, we may perform the detection with K different 
values r) v , v = 1...K, ranging from zero to a maximum value t\k = J] m ax equal to the nominal quantum 
efficiency of the detector. By denoting the probability of having m detected photons in the experiment with 
quantum efficiency -q = r) v by q™ = q™ we can rewrite Eqs. ® and (0 as 



- ^ A„(m,n)g n , (5) 



where v = 1, ..K and m = 0, .., M — 1. Let us now suppose that the £ n 's are negligible for n > N and 
that the r] u 's are known, then \fr} v Eq. © may be rewritten a finite sum over n from m to N. Overall, 
we obtain a finite linear system with K x M equations in the N unknowns {g n }. Unfortunately, the 
reconstruction of g n by matrix inversion cannot be used in practice since it would require an unreasonable 
number of experimental runs to assure the necessary precision HHH^l . This problem can be circumvented 
by considering Eqs. (J5]l as a statistical model for the parameters {g n } to be solved by maximum-likelihood 
(ML) estimation. The likelihood functional is given by the global probability of the sample i.e 

K M-l 

L=nn (^r- • (6) 

The ML estimates of {g n } are the values that maximizes L. In Eq. © n mu denotes the number of events 
"m detected photons" obtained with quantum efficiency r] v . The maximization of L under the conditions 
p n > O j g n = 1, can be well approximated by using the expectation-maximization (EM) algorithm 
EHH3 18], which leads to the iterative solution 



Q { n +1) =Q X n' I ;> ?^ A v [m 




-1 



m, n) v ^ . (7) 

„=lm=0 <l™[{Qn }] 

In Eq. (Q q$ denotes the n-th element of reconstructed statistics at the i-th step, 9jT[{^n^}] the theoretical 
probabilities as calculated from Eq. © at the z-th step, whereas f™ = n mu /n u represents the frequency of 
the event "m detected photons" with quantum efficiency r\ v , being n v the total number of runs performed 
with i] = rj u . In the following we assume equal n v , Vza 

Eq. (0 provides a solution once an initial distribution is chosen. In the following we will always consider 
an initial uniform distribution £^ = (1 + ./V) -1 Vn in the truncated Fock space n = 0, N. Results from 
Monte Carlo simulated experiments shows that any distribution with g$ ^ \/n performs equally well, 
i.e. the choice of the initial distribution does not alter the quality of reconstruction, though it may slightly 
affect the convergence properties of the algorithm. 



III. MONTE CARLO SIMULATED EXPERIMENTS 

In this Section we report the photon statistics for different states of a single-mode radiation field as 
obtained by ML reconstruction from Monte Carlo simulated photodetection performed with the same low 
quantum efficiency and different resolution threshold. We consider semiclassical states as well as highly 
nonclassical states and show results for different values of the maximal quantum efficiency 7? max - 




FIG. 1: Maximum-likelihood reconstruction of photon statistics using photodetectors with different counting capa- 
bility (different values of M) and maximum quantum efficiency r) max = 0.2. First line: photon distribution of a 
coherent state with (a^a) = 3; panel a) the reconstructed p n ; b): total absolute error eM; according to Eq. (|8j c): 
fidelity C?W according to Eq. l|9}. Results are reported for M = 1 (on/off detector), M = 2, M = 3 and M = 6. 
We use i-C = 30 different quantum efficiencies r) = r\ v , uniformly distributed in [q max /K, rj max \. The Hilbert space 
have been truncated at N = 30 and n v = 10 6 runs have been performed for each value t\ v . The last iteration is 
= n v . Other lines: same as first line for thermal state with (a^a) = 3, Fock state with n = 4; superposition of 
2 Fock states ( 1 3) + |5)). The inset of panel c) in the third line shows the last iterations in an expanded scale. 



In order to assess our results we introduce two figures of merit, to measure the reliability and the accuracy of 
the method respectively. Since the solution of the ML estimation is obtained iteratively, the most important 
aspect to keep under control is the convergence of the algorithm. A suitable parameter to evaluate the degree 
of convergence at the i-th iteration is the total absolute error 

K M 

^ =Y.H\^ - ^^}\\ ■ (8) 

Indeed, the total error measures the distance of the probabilities as calculated at the i-th iteration, 

from the actual experimental frequencies and thus, besides convergence, it quantifies how the estimated 
distribution reproduces the experimental data. The total distance is a decreasing function of the number of 
iterations. Its stationary value is proportional to the accuracy of the experimental frequencies {f™}. For 
finite data sample this value is of order for each value of r]^. As a measure of accuracy at the i-th 

step we adopt the so-called fidelity between probability distributions 

N— 



n=0 



In Ref. flill it was shown that using on/off detection (i.e. M = 1) a reliable reconstruction scheme may be 
obtained. In this paper, our aim is to check whether a higher value of M leads to some advantages, either in 
terms of accuracy or convergence. 

Fig. ^summarizes results for different states and values of M, assuming ?? max = 0.2 as the maximum 
value of the quantum efficiency. In order to compare the performances achievable by different M one has 
first to chose a criterion to stop the iterative algorithm. A natural criterion would be that of stopping the 
iteration when the total absolute error e^ 1 ' converges, i.e when the rate of its variation falls below a threshold 
and becomes negligible. This is also motivated by the behavior of eW versus the number of iterations (see 
central column in Fig. a rapid fall followed by a plateau. However, the most convenient value for the 
threshold unavoidably depends on the shape of the unknown distribution Qi n \ . Therefore, in the absence 
of any a priori information, this simple criterion should be supplemented by some additional recipes. We 
found numerically that a suitable criterion, valid for any value of M and, in the average, for any class of 
states, is the following: if the total absolute error appears to converge, then stop the algorithm at a number 
of iterations = n u equal to the number of measurements taken at each value of the quantum efficiency, 
otherwise stop the algorithm as soon as you see convergence. In cases when gW converges for a number of 
iteration < n v it is not convenient to stop the algorithm, since it may further increase the quality of the 
reconstruction, see for example G^> in the fourth line of Fig^ On other hands, if goes far beyond the 
condition of convergence of eW, then the algorithm may lose precision due to the noise of the data f™ (see 
inset in the second line of Fig. Q. The choice of the threshold at = n„ fits with other two independent 
facts. On one hand, we performed an extensive set of simulations and no considerable reduction of was 
observed before the convergence condition = n v . On the other hand, we applied the reconstruction 
algorithm to the exact on/off frequencies (not to simulated data) and no reduction of G^ 1 ' was observed for 
increasing #(2), thus confirming the data-noise origin of precision loss. 

Let us now illustrate our results about the accuracy of the reconstructions. As mentioned above we 
performed simulated experiments, on semiclassical states like coherent and thermal ones as well as highly 
nonclassical states such Fock states, by using three values for the maximum quantum efficiency r/ max = 
0.2, 0.5, 0.8. Simulations have been performed by using n v = 10 4 data at each value of the quantum 
efficiency. This value of n„ is relatively small and well within the realm of quantum optical experiment. In 
a real experiment performed under the same conditions, n v of this order would allow on-line reconstruction 
of the photon statistics. Simulations performed using different values of the parameters lead to similar 
conclusion. 

In Fig. |2 we show the fidelity of reconstruction £?M a t the last iterations, G^ L \ for different input 
signals and different resolution thresholds as a function of the average number of photons of the signal. As 
it is apparent from the plots, for low quantum efficiency and semiclassical states on/off detectors (M = 1) 
are enough to achieve a good reconstruction. On the other hand, the accuracy for nonclassical states largely 
improves for M > 2. For higher values of the quantum efficiency (the middle and the right plots) on/off 
detectors becomes sufficient for a good reconstruction also for nonclassical states having single-peaked 
distributions. Notice that a higher resolution, however with a low value of rj max , does not guarantee a good 
reconstruction. Also notice that the results reported in Fig. [2] has been obtained by stop the algorithm 
according to the prescriptions mentioned above. We therefore do not expect that the fidelity is optimal for 
any state. Indeed the fidelity G^ L > in panel C5) and C8) of Fig. [2] slightly decreases for high number of 
photons, though remaining close to unit value. The decreases in panels T2), T5) and T8) of Fig. [2]for greater 
mean values is due to the choice of a small dimension of the truncated Hilbert space: by increasing N the 
trend disappear. The low values of G^ %L ' for the superposition states are due to the noise in f™, also for the 
high quantum efficiencies. Indeed, by increasing n v the reconstructions greatly improves; besides, for high 
quantum efficiency the detectors perform good reconstructions also for M = 0, see the inset of panel S2) 
in Fig. 12 The error bars in Fig. |2]are standard deviations of G^ L ^ as calculated from 40 different Monte 
Carlo runs. They may appear large, but this is due to the high value of G^ L \ which in turn determines the 
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FIG. 2: Fidelity parameter at the last iteration, G^ lL \ as a function of (a^a) for 4 different statistics of coherent, 
panels labeled C); thermal, label T); Fock state, label F); superposition of 2 Fock states, -^(| (a^a) — 1) + 1 (at a) + 1)), 
label S). The reconstructions have been performed using r) max = 0.2, 0.5 and 0.8, panels labeled with 2), 5) and 8) 
respectively; n v = 10 4 number of runs at each for every r\ v . According to the rule given in the text, the last iteration 
performed is = n„ except for coherent state with M = 1 and (a^a) = 1...4 for which = 10 5 , and for 

the superpositions state. Inset of S2) show G^ L ^ for ^(|5) + |7)), M=0, quantum efficiency 0.2 and 0.8 in function 
of n v . Every point of the curves is the mean of 40 simulations, the error bars represent the standard deviations. The 
other parameters are the same as in fig[D 



scale of the plot. Look at panels S)' for comparison. 



IV. EXPERIMENTAL DATA 

In order to confirm the Monte Carlo results for single-peaked distributions we have performed the recon- 
struction of the photon statistics of a coherent signal obtained from a Nd:YLF laser. The experimental data 
have been recorded with a hybrid photo detector, Hamamatsu H8236-40, placed on the second harmonics 
(523.5 nm) of a cw mode-locked Nd:YLF laser regeneratively amplified at a repetition rate of 5 kHz (High 
Q Laser). 

The frequencies f™, until m = 3, have been extrapolated from the response of the detector with the 
following procedure. At first a Gaussian best-fit of the response peaks at 0,..,4 photoelectrons has been 
performed 10], then the the frequencies has been obtained by choosing three thresholds, whose optimal 
choice turns out to be the mean point between photoelectron peaks II 911 . After the frequencies f™ have 
been obtained we used the algorithm of the previous Section. The results of the reconstruction are shown 
in Fig. |3j panel a). The reconstructed distribution at the last iteration Qn L ^ have been then compared with a 
Poissonian best-fit. As it is apparent from Fig. |3] the agreement between Qn L ^ and the best-fit is very good. 



In addition, we notice that the shape of C?W and are very similar to those coming from simulation and 

(i) 

that the g n obtained by the experimental data with M=l is close to the one obtained with M=2 and M=3, 
as predicted by simulations. 
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FIG. 3: Experimental reconstruction of the photon number distribution of a coherent state, a): reconstructed p n and 
coherent best fit, (a^a) = 2.98; b): total absolute error e-W; c): fidelity calculated a posteriori. The reconstruction 
is performed with a hybrid photo detector operated by taking M = 1 (on/off detector), M = 2 and M = 3; K = 100 
different quantum efficiencies 77 = ?/„ distributed in [0, T) max \, rj max = 0.4; the Hilbert space is truncated at N = 30; 
n u = 10 4 number of runs have been performed for each 77; The algorithm is stopped at iteration = n v . 



V. CONCLUSION 

We have compared the reconstruction of the photon statistics as obtained applying ML algorithm to data 
coming from detectors with different counting capability and different quantum efficiencies. We found that 
using ML methods, detectors with high quantum efficiency does not need to have high counting capability, 
since on/off detection assisted by ML methods already provides good state reconstructions. On the other 
hand, a small quantum efficiency makes the counting capability a crucial parameter. Overall, our results 
indicate that development of future photodetectors may be focused on increasing the quantum efficiency 
rather than the counting capability. 
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